import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt


def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y


T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T - 1:
    if t - 2 < 0:
        t=t+1
        continue
    pi_star = norm.rvs(loc=pi[t - 2], scale=sigma, random_state=None)
    print (11,pi_star)
    # 评价指标是mean3 sigma2的正态分布概率
    a = norm_dist_prob(pi_star) / norm_dist_prob(pi[t - 2])

    alpha = min(1, a)
    u = random.uniform(0, 1)
    if u < alpha:
        pi[t - 1] = pi_star
    else:
        pi[t - 1] = pi[t - 2]
    t = t + 1
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
plt.hist(pi, 50, facecolor='red', alpha=0.7)
plt.show()
